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Discontinuous Galerkin (DG) methods coupled to WENO algorithms allow high order convergence 
for smooth problems and for the simulation of discontinuities and shocks. In this work, we investi¬ 
gate WENO-DG algorithms in the context of numerical general relativity, in particular for general 
relativistic hydrodynamics. We implement the standard WENO method at different orders, a com¬ 
pact (simple) WENO scheme, as well as an alternative subcell evolution algorithm. To evaluate the 
performance of the different numerical schemes, we study non-relativistic, special relativistic, and 
general relativistic testbeds. We present the first three-dimensional simulations of general relativis¬ 
tic hydrodynamics, albeit for a fixed spacetime background, within the framework of WENO-DG 
methods. The most important testbed is a single TOV-star in three dimensions, showing that long 
term stable simulations of single isolated neutron stars can be obtained with WENO-DG methods. 
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I. INTRODUCTION 

Over the last decade simulations in numerical general 
relativity have seen a tremendous improvement in accu¬ 
racy and stability and have become an important tool 
for the study of high energy and strong gravitational 
field effects. To date numerical simulations are the only 
possibility to investigate complex astrophysical scenarios 
as e.g. stellar collapse [T] and coalescing binary neutron 
stars [2] ■ Although numerical simulations are in principle 
not restricted by approximations beyond the numerical 
approximation, they are limited by the finite accuracy of 
the particular discretization method. Among the differ¬ 
ent methods to solve partial differential equations like 
those of general relativity, the discontinuous Galerkin 
(DG) method has emerged in recent years as a partic¬ 
ularly successful general purpose paradigm I5H5]. It can 
be argued that the DG method, more explicitly the DG 
finite element method or DG spectral element method, 
subsumes and combines several of the key advantages 
of traditional finite element and finite volume methods, 
e.g. m. In particular, the discontinuous Galerkin method 
works with element-local stencils, which is a great advan¬ 
tage for parallelization and the construction of compli¬ 
cated grids. Furthermore, DG methods offer easy access 
to hp-adaptivity [B] , where both the size of the computa¬ 
tional elements (or cells) and the order of the polynomial 
approximation within each element can be adapted to the 
problem. For smooth solutions, DG methods approach 
the optimal order of exponential convergence of pseu- 
dospectral methods on multiple patches. In fact, certain 
DG methods are equivalent to pseudospectral methods 
with a specific penalty method for the patch boundaries 
[7]. For non-smooth solutions, low order elements have 
been combined with various HRSG schemes, for example 
in the form of WENO DG methods [5]. 


In this work we consider the application of DG meth¬ 
ods to simulations in numerical general relativity cou¬ 
pled to general relativistic hydrodynamics (GRHD). Gon- 
cretely, the goal is to compute the numerical evolution 
of spacetimes containing neutron stars. The governing 
differential equations are the time-dependent, non-linear 
Einstein field equations for the spacetime geometry cou¬ 
pled to a relativistic fluid model. Most relativistic hydro¬ 
dynamics simulations are based on the “Valencia formu¬ 
lation” , in which the matter field evolution is given in a 
conservative form [^. 

Among the numerous numerical studies carried out in 
the field, most have been performed using finite difference 
(FD) and finite volume (FV) methods, with significant 
success. For the geometry (including black holes), high- 
order finite differencing is the rule, often 4th to 8th order 
finite differences in space for structured adaptive mesh 
refinement (AMR), e.g. [TUIITT] . The matter part allows 
the formation of strong relativistic shocks, and a vari¬ 
ety of finite volume (or finite difference) HRSG schemes 
have been developed mm- For smooth solutions, pseu¬ 
dospectral methods have been very successful [T51116) . 
Recently, a convergence order of ~ 3 was observed for 
high order matter formulations in [IIHIS]- 

DG methods for numerical relativity offer the usual list 
of attractive features. In particular, one goal would be to 
combine high-order, smooth regions with lower-order re¬ 
gions containing shocks. Compared to AMR with large, 
overlapping finite difference stencils, the DG spectral el¬ 
ement method is more easily and more efficiently paral- 
lelizable, while still allowing high-order approximations. 
However, there remain several open issues with regard 
to DG methods in numerical relativity. Some issues are 
known but unresolved, some have simply not been inves¬ 
tigated yet. 

The evolution equations of numerical relativity are a 
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coupled system for the geometry (the metric variables) 
and the matter variables. While the matter equations 
are naturally given in a flux form [^, this is not the case 
for the geometry. Since a typical DG method starts with 
a flux-balance law, it is in principle straightforward to 
design a method for the matter part. On the other hand, 
for the geometric part one should either recast the equa¬ 
tions in a hyperbolic flux form, or suggest less standard 
methods. 

There have been essentially only three major efforts to 
employ DG methods for general relativity and/or GRHD. 
In [20] , Zumbusch gives the first and so far only example 
for a complete DG method for the (3 -I- l)-dimensional 
(short 3D) Einstein equations in vacuum. Discussed is 
a space-time DG scheme mostly in the context of lin¬ 
earized equations and in a specific gauge, but the scheme 
also handles non-linearities. So far there has not been an 
astrophysics application, say involving black holes or neu¬ 
tron stars. In ED, Brown et al. discuss a DG method for 
the so-called BSSN formulation of the vacuum equations, 
mostly with ID examples. The BSSN equations are not 
in flux-form, but the various non-linearities and second 
derivatives are successfully dealt with on a case by case 
basis. And in [22], Radice and Rezzolla give a general 
discussion of the DG method for the matter equations 
in the standard flux form, without including the equa¬ 
tion for the geometry. They present a working ID im¬ 
plementation for general relativistic matter in spherical 
symmetry. In addition, there has been work on special 
relativistic hydrodynamics (SRHD). Zhao and Tang [25] 
were the first to apply the WENO-DG method of [5] to a 
variety of ID and 2D test cases in SRHD, and the method 
turned out to be robust and reliable in capturing shocks. 

The concrete target of the present work is to model a 
single stationary neutron star (a TOV star [211 |2S|) in 
3D, although we perform a variety of tests in ID and 
2D as well. The TOV star is computed in the Gowling 
approximation, which simplifies the problem by assuming 
that the geometry may be curved but does not depend on 
time, which in turn is compatible with the stationarity 
of the TOV star. The numerical evolution of the matter 
variables for fixed metric is a standard approach that 
still allows to test key features of the hydrodynamics, 
including the treatment of the non-differentiable density 
at the surface of the star. We leave the coupling to a 
dynamic geometry to future work. 

In preparation for the simulations in full, 3D GRHD, 
we test the Runge-Kutta DG (RKDG) method coupled 
to a variety of WENO reconstructions for the equations of 
general relativistic hydrodynamics jS] . We reproduce the 
non-relativistic standard results [HI US] as well as some of 
the special relativistic test cases of [23] for a third and 
a fifth order method, WEN03 and WEN05. We extend 
[25] by also considering WENO-Z m and the simple 
WENO limiters of [21]. Finally, we present the first ap¬ 
plication of RKDG WENO methods to a 3D TOV star in 
the Gowling approximation. The numerical experiments 
are implemented in the new bamps code m for spec¬ 


tral element methods. We import some methods from an 
existing full-featured finite difference AMR code for 3D 
GRHD, BAM [Tni[2HH3I]- 

When researching the available HRSG methods for 
DG, there is one issue related to shock resolution and 
efhciency that is well-known but that does not always 
appear to receive the attention it deserves. For ED or 
FV methods with HRSG, shocks are resolved within a 
few cell widths, which means within a few grid points. 
For DG methods, the standard approach is to employ 
WENO reconstruction based on cell averages [8]. In such 
WENO-DG methods, shocks are again resolved within a 
few cell widths, but each cell now contains p points (for 
polynomials of order p— 1). The WEN03 stencil involves 
3 cells and 3p points, and the WEN05 stencil involves 
5 cells and 5p points. Effectively, the high resolution 
within each cell (the “subcell resolution”) is lost if only 
the cell averages are used. For practical implementations 
a rough estimate is therefore that such WENO-DG meth¬ 
ods could require about 27 or 125 times more resources 
for shock resolution in 3D than comparable FD or FV 
methods (these factors vary with the actual implementa¬ 
tion). 

For the evaluation of DG methods for GRHD it mat¬ 
ters whether such methods are competitive to existing 
FD/FV methods in terms of efficiency. Hence we consider 
the following measures aimed at handling the compara¬ 
tively low efficiency of cell-averaged WENO-DG meth¬ 
ods. A common strategy is to limit the application of 
the WENO scheme to only those cells that need it, and 
there has been quite some work on so-called “troubled 
cell indicators” [52] . 

A recent development are the so-called “simple” 
WENO methods of [26], which effectively construct a 
compact stencil for high-order WENO methods. For ex¬ 
ample, the fifth-order WENO method is constructed from 
only 3 instead of 5 cells, using the high-order informa¬ 
tion from the nearest neighbor cells to obtain fifth order. 
This leads to significant savings, but the method has not 
been widely tested yet. We include the compact/simple 
WENO method in our tests and report on some differ¬ 
ences to the standard WENO method, in particular in 
3D. 

Another important development is a hybrid approach 
[55H55] . which replaces troubled cells by an equidistant 
subgrid and applies FV shock capturing on these grids. 
This approach maintains the subcell resolution of FV 
methods, but increases the complexity of the implemen¬ 
tation since two types of grids and special grid transfer 
operators are required. In our case the method is appeal¬ 
ing because a full-featured FD implementation is already 
available [2H1 [12]. If successful, the strategy would be 
to construct a high-order DG method for regions where 
the solution is smooth, but to rely on established FD 
methods near shocks. 

HRSG for DG comes at a cost since WENO-DG as well 
as the hybrid FD-DG method break the cell-locality of 
the basic DG method. We consider both methods here 
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to gain some insight into their relative merit. 

In Sec. [nj we introduce the DG method for 3D flux- 
balance laws, specify the equations of relativistic hydro¬ 
dynamics, and discuss the WENO-DG and FD-DG meth¬ 
ods. In Sec. m we summarize the numerical implemen¬ 
tation. As basic tests we consider the advection equation 
and the Burgers equation in ID in Sec. IV while ID and 
2D tests for SRHD are presented in Sec. |V| The main re¬ 
sults concern the evolution of a TOV star in Sec. IVD We 
conclude in Sec. m For completeness, we collect some 
relevant details of the basic ID DG method in App. fA} 

Throughout the article dimensionless units are used, 
i.e. we set c = G = Mq = 1. We denote spacetime indices 
by a, 6 ,... and indices over space dimensions by i,j,.... 


II. METHODS 

A. Discontinuous Galerkin method 

The hydrodynamical equations governing the time evo¬ 
lution of the matter fields can be cast as a non-linear con¬ 
servation law for a vector of variables, u(a;, t), depending 
on time t and position x G K^. The conservation law is 
given by 


5tu-|-air(u) = S , (1) 

with the sources S and the fluxes f®. We summarize some 
of the relevant aspects of the DG method for conservation 
laws of scalar function on K in Appendix]^ while simply 
stating the key equations for vector-valued functions on 
K" here, cmp. [A 11 ^. 

We consider a partition of M" into cells Ij, x G Ij , and 
define the finite dimensional approximation space 


at the boundary, u_ and u_|_, and reproduce the original 
flux if u is continuous. A simple example of a numeri¬ 
cal flux with this property is the local Lax-Friedrich 
(LLF) flux 

f**(u_,u+)ni = i [r(u_)ni -H r(u+)ni - A (u+ - u_)] 

(4) 

where A denotes the maximum absolute eigenvalue of 
the Jacobian d{Pni)/du. We use an LLF algorithm 
throughout this article. Writing out the numerical so¬ 
lution u(x, t) as an element of explicitly, 


N 


u\iAx,t) = Vufe(t)u'"(x) 


A;=0 


(5) 


and a basis of P^(Ij), allows to recast ([^ as an al¬ 
gebraic equation for the unknown time derivatives 9tUfe. 
To evolve these coefficients in time, we use an explicit 
fourth order Runge-Kutta method. More details on the 


actual implementation are given in Sec. Ill 


B. Relativistic hydrodynamics 

Although we are working in cowling approximation, 
i.e. keeping the metric fixed, the matter fields are evolved 
dynamically on a curved spacetime background. We want 
to recast briefly the important equations and methods 
necessary to solve the general relativistic hydrodynamical 
equations; special relativity can be easily obtained by 
choosing flat spacetime. 


1. 3+1-decomposition 


F^:={u:u(x)|z, GP^(J,)} (2) 

with P^(/_,) denoting the finite dimensional space of 
polynomials on Ij of degree at most N. As in most 
standard applications, we set the polynomial order N 
as a constant over the whole partition. To deduce a DG 
scheme from Eq. Q, we want to find a function u„(x) 
for which the weak form 

f vdtUn dV + f r(u„)vni dS” 

Jlj Jdl, 

- [ P{un)diV dV= [ Sv dF (3) 

Jij Jij 

holds for all v G . For simplicity, we denote the ap¬ 
proximate/numerical solution Un(x) as u(x) in the fol¬ 
lowing. An important advantage of the DG-scheme is 
that V does not need to be continuous at the cell bound¬ 
aries. Therefore, no unambiguous definition of the fluxes 
at cell boundaries entering Eq. (|^ exists. To overcome 
this issue, we introduce the numerical fluxes f**(u_, u+), 
which depend on the inside/outside cell limited value of u 


Although we assume the spacetime to be fixed, we have 
to recast it in a suitable form for dynamical evolutions. 
This can be done with the help of a 3 -I- 1 decomposi¬ 
tion [371 [3H] (see |53MT] for textbook introductions) in 
which the four-dimensional spacetime metric is rewritten 
as 

ds^ = —a^dt^ + 7ij(dx* -I- /3Mf)(dx'^ -|- P^dt), (6) 

where a is the lapse function, /3® the shift vector, and 
the spatial metric. In case of a flat spacetime a = 1, /3* = 
0, 7 ij = 6ij employing Cartesian coordinates. Einstein’s 
field equations split into two sets, the constraint equa¬ 
tions and the evolution equations. For our single neutron 
stars tests, we recast the Tolman-Oppenheimer-Volkoff 
(TOV)-equation [331|33 in 3-1-1-form and solve it to ob¬ 
tain an ordinary differential equations. In addition to the 
3-|-1-split we perform a conformal transformation of the 
spatial metric, 

(7) 

where tp is the conformal factor and 7 ^ the conformally 
related metric. 
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2. Hydrodynamic equations 

According to Eq. Q we denote the state vector collect¬ 
ing the conserved variables as u, while f*(u) are hydro- 
dynamical fluxes, and S the source terms. The fluxes and 
the sources depend in general on the metric and matter 
flelds. The conserved variables are u = Sk, t), 

and denote respectively the rest-mass density (D), the 
momentum density (Sk), and an internal energy (r) mea¬ 
sured by the Eulerian observer given by the particular 
spacetime foliation. 7 = det "fij is the determinant of the 
spatial three-metric. The conserved variables u can be re¬ 
constructed from the primitive variables w = (p, u*, e,p), 
i.e. rest-mass density, 3-velocity measured by the Eule¬ 
rian observer, internal energy and pressure of the fluid 


by the following equations: 


D = Wp, 

( 8 a) 

Sk = W^phvk, 

( 8 b) 

T = {W^ph — p) — D, 

( 8 c) 


where W is the Lorentz factor, W = l/Vl — ViV'^ and h 
is the specific enthalpy h = 1 + e + p/p. 

To close the system an equation of state (EOS) p = 
P{p, e) is needed. In this work, we use a simple polytropic 

P{p) = Kp^ (9) 

or an ideal gas EOS of the form 

Pip,e) = {T-l)pe. (10) 

The particular implementation of the hydrodynamical 
equations follows [gEH]. 

However, due to the special choice of the background 
metric, the flux and source terms simplify dramatically 
by setting jij = Sij and /3* = 0 in all our examples. 


3. Primitive recovery and atmosphere treatment 


We evolve the conservative variables u by constructing 
the fluxes and source terms for every time slice. While 
f* and S both contain the primitive variables w we have 
to recover those from the conservatives. The inverse re¬ 


lations of ( 8 a)-( 8 cl are given by: 


P = 


e = 


D 

W’ 

5 * 

T + D +p^ 

\/{t -\-p-\- DY — 5'^ — Wp — D 


( 11 ) 

( 12 ) 

(13) 


with W = (t + p + D)/ + P + Py — S'^ and = 

SiS^. To make use of (11)-(13|, we have to determine 
the pressure p. 

The explicit primitive reconstruction goes as follows. 
First, we try to recover the primitive variables for the 


full equation of state including thermal components p = 
P{p,e). For this reason a Newton-Raphson method is 
employed to compute the pressure p. If the method does 
not converge to the desired accuracy a cold equation of 
state p = p{p) is used and we try to find with a Newton- 
Raphson method the density p. 

As in most general relativistic hydrodynamic codes, we 
have to include an artificial atmosphere to solve the prob¬ 
lem of fluid-vacuum interfaces. This allows long term 
stable and robust numerical simulations . The at¬ 

mosphere Patm is computed according to 

Patm = /atm • max[p(t = 0)]. (14) 

Whenever a point falls below the atmosphere threshold 
Pthr = /thr ’ Patm during the evolution or the primitive 
reconstruction, it is set to the atmosphere value. 


C. WENO reconstruction methods 


As a next step, we explain how to avoid oscillations and 
unphysical behavior caused by the Gibbs phenomenon. 
For this purpose we locate discontinuities and oscillations 
with the troubled cell indicator described in Sec. IHCTI 
and apply a WENO limiter reconstruction [HIES!. We 
introduce three different WENO reconstruction meth¬ 


ods, the standard WENO approach (Sec. IIC 2), the 
simple WENO algorithm [26] based on compact stencils 
(Sec. lie3), a nd a W ENO algorithm based on a subcell 
evolution (Sec. IID). 


1. Troubled cell indication 

Given the coefficients of the numerical solution Up(t) 
at time t, we can calculate the average of the polynomial 
u{x,t) over the grid patches Ij = [oj, 6 j]: 

Uj := ^ / Up(/>^(a:) dx=]^ J ^ ■ (15) 


We further denote the boundary values of u as 


Uj :=u(aj). 

:=u(6j) 

(16) 

and define the four differences: 



1 

'll 

1 

c 

1 

■= < - "j 

(17) 

A_u:= Uj-Uj_i, 

A+u := Uj+i - Uj 

(18) 


We also introduce the minmod function 
minmod(a;i, a: 2 ,..., Xn) = 

I s ■ mini<j<„ \xj\ if sign(a:i) = ... = sign(x„) =: s 
1 0 otherwise 

(19) 
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FIG. 1: The WENO-5 (w = 2) methodology applied in a smooth case (left figure) and a shock case (right figure). The values 
in the interval x £ Ij = [—1,1] are to be reconstructed from the five grid patch averages Uj- 2 ,Uj-i,Uj,Uj+i,Uj+ 2 - The three 
stencils Si, 82,83 are created as a clustering of three grid patches each with the corresponding approximating polynomial 
pi{x),p 2 {x),P 3 {x). Another higher order polynomial q{x) can be found from employing all hve averages. Following the strategy 
as described in |II C| the smoothness indicators j3i are calculated for each stencil. A large j3i indicates non-smoothness of the 
corresponding polynomial Pi, which leads to a minor contribution of the stencil Si for the reconstruction. In the shock case, 
the reconstructed point values (empty black circles) lie very close to the smoothest polynomial ps, whereas in the smooth case 
all three approximating polynomials are taken into account almost equally, so that the reconstruction is very close to the 5th 
order polynomial q (filled gray circles). 


and the modified TVB minmod function 


in a:-direction 


minmodTVB(a;i,a; 2 , ...,a;„) = 

{ oi if |oi| < M (maxj Axj)"^ 

minmod(xi,a; 2 , ...,x„) otherwise 


( 20 ) 


' AyAz 

u-*- ■=^ 
■ AyAz 


nbk rbi 

ak ^ ai 
pbk pbi 


\ipqr(t^{aj)(j)^{y)(j/'{z) Aydz, 


Upgr^{bj)(l)‘^{y)(j)^{z) dydz. 


( 22 ) 


Our troubled cell indicator marks a grid patch as trou¬ 
bled, if 


minmodTVB ((mj , A_u^, A+u'^^ ^ {uj or 
minmodTVB ^ (“/)* (21) 


for at least one component k. This is exemplary for a 
situation, in which a component of u is not monotonous 
(because the arguments of minmod differ in sign) or its 
gradient inside a patch is larger, than that of the neigh¬ 
boring patches (shock inside the cell). 


In the case of multiple dimensions, we perform the 
ID troubled cell indication in every direction. A cell is 
marked as troubled, if at least one of these indications re¬ 
sults in a troubled state. To apply the ID algorithm, the 
boundary values used in (16) have to be modified, since 
the cell boundaries are not longer single points, but lines 
or surfaces. Therefore, we redefine u)" by the boundary 
averages, i.e. for a 3D cell I = [aj,bj] x [ak,bk] x [ai,bi] 


2. Standard WENO reconstruction 

In a standard WENO method of order 2'u; -|- 1, 
we construct w -I- 1 stencils Si around Ij, each 
as an aggregation of re -I- 1 grid patches: Si = 

S: 't A tu. In Fig. ^ this 
partitioning is shown for w = 2. For each stencil, we con¬ 
struct a w-th order polynomial Pi, which has the same 
average as the numerical solution u over each grid patch 
in the stencil. That means solving the system 

Uk = ^ / Pi (a:) dx, for all 4 S Si (23) 

Ax Ji^ 

for the w + 1 coefficients of each component of Pi. Simi¬ 
larly, we construct a 2w-th order polynomial q fulfilling 

Uk = ^ / q(x) dx, for all 4 G S, (24) 

Aa: 

with S := UiSi being the large stencil over all 2w-|-l grid 
patches. The fundamental concept is to approximate the 
solution in [—1,1] as a linear combination of the Pi, which 
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should give the same result as the higher order approxi¬ 
mation q in smooth regions. This condition defines the 
linear (or ideal) weights 7 ^ satisfying 


W-\-l 


q(a;) = X! 


(25) 


i=l 


We emphasize that the 7 ^ depend on the point x where 
the approximation should hold. It is remarkable that al¬ 
though both sides of Eq. (251 depend intrinsically on the 


2 w + 1 averages and the system is overdetermined 
(only w -f 1 variables), we could always find an exact 
solution for ( |25[ ) in our tests. In regions where the so¬ 
lution is not smooth, the weights should be chosen such 
that the smoothest polynomial of Pi is preferred. For 
this purpose, we use a smoothness indicator as suggested 




1^1 




Because ,3^ is large for non-smooth pi, the weights are 
chosen indirect proportional to /3j. We use either the 
standard WENO choice 


U}i{x) = 




\2 ’ 


(10-6+ /3J^ 

or the improved WENO-Z version [27] for w = 2, 

1 / 33 -/3i I 


d>j(a;) = 7i(a;) 1 -f 


l3i + 10-40 


and normalize the results: 


Cj.ix) 




(27) 


(28) 


(29) 


where u}i(x) are the final reconstruction weights. The 
reconstructed solution is then given by: 


W-\-l 


uWENO^^.) ^ U}^{x)pi{ 


(30) 


To generalize the presented reconstruction mechanism 
to 2D and 3D, we use the procedure described in [23]. For 
simplicity, we assume a rectilinear 2D grid structure with 
iV -I- 1 grid points per cell and direction. To reduce 
the full reconstruction of the cell Ijk to the ID case, we 
decouple the different directions. First we perform 2w + l 
ID WENO reconstructions in the x direction with input 
data 




w,k’ 


’ 4^j+w,k|' 


k — w<k<k + w 

(31) 

to reconstruct the -|- 1 averages per cell: 

uj’- := f u{^p,y,t) dy, k-w<k<k + w, (32) 

1 < p < m -|- 1 


Then, we can apply a second ID WENO reconstruction 
based on the ID averages in y direction with the input 
data 


{upk-w.U? 


k—w-f 1 ’ 


.Uyk+w}> l<p<m + l 


(33) 


to get the 2D reconstructed values inside the cell Ijk- 


WENO 
■^jk 


iCpAq), l<p,q<m + l 


(34) 


3. Simple WENO reconstruction 


To reconstruct the polynomial with a standard WENO 
method, the cell averages of many neighboring cells are 
needed. This leads to large computational costs and an 
undesirable smoothing of the solution. In [26], it is dis¬ 
cussed, that this standard procedure is not necessary in a 
DO method, since the neighboring cells yield more infor¬ 
mation than only a cell average value. This idea implies 
the simple WENO reconstruction, in which the standard 
WENO methodology is applied to the cell polynomial 
and the neighboring cell polynomials by redefining 

p*(x) := u|/^_^.(x) -f (35) 

^ / (u|T (a^)-u|7^.+,(a;))(ia;, i =-1,0,1. 

The integral values cause a shift of the polynomials, so 
that they all have the same average value in cell Ij and 
the cell average is conserved during reconstruction. The 
corresponding expressions in (26) and ( j^ have to be 
substituted. Furthermore, we only use the next neigh¬ 
bors, which is setting w = 1 {3 cell stencil) in all WENO 
formulas. Since with the new ansatz ( ]35| ) every linear 
combination of the Pi(a:) is a higher order approxima¬ 
tion, there is no need to find special ideal weights as in the 
standard WENO method. Instead, we can freely choose 
the weights for all involved cells. In our tests, we choose 
7_i =71 = 1-10-6, 7 q = 1 — 271 = 1 — 2-10-6 for smooth 
setups and 7 _i = 71 = 1 - 10 -^, 79 = 1 —271 = 1 - 2 - 10-6 
for problems with discontinuities. 


D. Subcell evolution method 

Finally, we consider a hybrid FD-DG method moti¬ 
vated by [33] where shock capturing is performed on a 
subgrid of equidistant grid points. The method of [33] 
is based on subgrids, an a posteriori troubled cell indica¬ 
tor, and a locally implicit time integrator. We decided to 
investigate the subgrid method separately without these 
other features, so we cannot compare directly to [33] . 
There are open questions regarding the stability and ac¬ 
curacy of the subgrid method, in particular when used 
without the locally implicit time integration method (but 
see also [3S1|35]). For our method we use the same trou¬ 


bled cell indicator as introduced in IIC If a cell Ij has 
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been marked as troubled, we subdivide this cell in 2N +1 
equidistant subcells Jk containing a single point yk each 
(where N is the polynomial order) and compute the value 
of the approximating polynomial on the individual sub¬ 
cells points: 

vk = Up(()^(yfc), for all G. (36) 


successfully used to study the gravitational wave col¬ 
lapse and it allows long-term simulations of single black 
hole spacetimes with excision techniques. The program 
exhibits a hybrid p-thread/MPI parallelization strategy 
and shows almost ideal scaling for up to several thou¬ 
sands of computing cores in vacuum simulations, see [14] 
for more details. 


This map Up i—^ Vk can be done with the subcell projec¬ 
tion operator V. The back projection 'P~^ is non-trivial, 
because the problem of finding a polynomial of order N to 
satisty the given 2N +1 equations (36) is overdetermined. 


Performing a least-squares fit of a 7V-th order polynomial 
for the 2N + 1 points turns out to be a good choice for a 
back projection. In our tests, we found the correspond¬ 
ing matrices for V and V~^ to be pseudoinverse. This 
is easy to verify, because whenever Vk originate from an 
exact fV-th order polynomial, a least-squares fit 7^“^ will 
give the exact polynomial coefficients, so 'P~^'P = 1 (but 
not neccessarily 'PV~^ = !)■ It is important to notice 
that contrary to |33j . we use a projection matrix based 
on the point values in the subcells, not the averages Vk. 
This is necessary, since we want to employ a FD code on 
the subcells instead of a FV method, leading to a viola¬ 
tion of conservation laws (e.g. of the rest-mass), when a 
projection from topcell to subcells, or vice versa is done. 
However, in our tests we found this defects decaying with 
order TV -|- 1, when we raise the grid resolution. We im¬ 
port all necessary routines of the BAM code [101 EH ES] ■ 
In [551 05] this scheme is explained in detail. Further 
improvements allow to obtain high order convergence in 
smooth regions will be presented in [46]. The general 
idea is to discretize Eq. Q as 


9tUk = 


2N + 1 
Ax 


k-. 


- F 


k-l-l 


(37) 


with Ax being the cell grid spacing, 2 n+i subcell grid 
spacing and Fk+i the numerical flux at the boundary 
between subcells Jk and Jk+i- 

The subcell interface values of the fluxes fkG 1 are com¬ 
puted with the LLF scheme. The necessary right and left 
states for the interface flux calculation are provided by 
a WENOZ [^03] reconstruction from the given subcell 
values. Having evaluated the RHS of Eq. (37), we use an 


explicit fourth-order Runge-Kutta method for the time 
step. After each Runge-Kutta substep the new subcell 
values are back projected to the DG-grid by . If the 
cell stays troubled in the next time step, the next evo¬ 
lution step is based on the subcell results without using 
the back-projected results. 


In this work we extend the bamips code by imple¬ 
menting (i) discontinuous Galerkin methods, (ii) a gen¬ 
eral relativistic hydrodynamics scheme for fixed back¬ 
ground metrics, (iii) a simple high resolution shock cap¬ 
turing (HRSG) scheme as in [23], (iv) a subcell-HRSG 
scheme [33]. This work is the first step towards a more 
general infrastructure for the simulation of compact bi¬ 
nary systems where matter is present. 

Although bamps allows grid structures known as 
’’cubed spheres” [47], we restrict ourselves to simple 
Cartesian boxes. However, a generalization could be 
achieved easily. For the actual implementation of Eq. (§, 
we map each Ij to a reference box [—1,1]^ and define 
N + 1 Legendre-Gauss-Lobatto points G [—1,1] for 
each direction. Given these points, we choose the basis 
yk of P^([—1,1]3) to be the product of the corresponding 
Lagrange interpolating polynomials each applied to one 
component of ^ 


= 


with 


- n IfI 

j^p 


(38) 


(39) 


i.e. we use a nodal DG formulation. The chosen basis 
allows us to use = 5^'^ and simplifies the computa¬ 
tion of the coefficients \ipqr = u(^p, ^’') (interpolation 

condition). In contrast to the modal DG formulation, the 
flux and source coefficients are then easily determined by 
pointwise evaluations = f®(upgr)) Spgr = S(upgr-)- 
Defining the mass matrix 


M^b 





(40) 


and the stiffness matrix 


III. NUMERICAL IMPLEMENTATION 

Throughout this article we employ the bamps code [14]. 
It is based on the method-of-lines with a pseudospectral 
decomposition in the spatial part and an explicit fourth- 
order Runge-Kutta for the time stepping. It has been 


5 '“'’ = J (41) 

we seperate analytic expressions and numerical variables 







in Eq. ([^ to gain the semidiscrete scheme: 


^t^pqr 


^bqr 


^ oa.bf2 

Ay V 9“ 

- i,r) - 

- e, 1) - e, -1) 


Spqr ■ 


(42) 


Due to the choice of collocation points, the mass and stiff¬ 
ness matrix can be determined using Legendre-Gauss- 
Lobatto integration with the corresponding weights 


gab ^ 


(43) 

(44) 


Notice, that Eq. (43) is just an approximation, while (44) 
is exact, since the N + 1-point Legendre-Gauss-Lobatto 
integration is exact for polynomials of order 2N— 1. This 
approximation simplifies the scheme and brings M in a 
diagonal form. Furthermore, it is equal to a modal filter, 
which decreases the highest mode by a factor N/ (27V -|- 

1 ) m- 

For the standard WENO implementation, we recast 
the crucial equations in matrix form, where all matrices 
can be precomputed from the geometry before evolution. 
During the actual simulation (i) the smoothness indica¬ 
tors are calculated from the cell averages as a quadratic 
form Pi = Qf‘Uij^kUi+i', (h) the weights are determined 
by ( p^ ; (iii) the value Pi{^‘^) of the approximating poly¬ 
nomial of stencil i at the collocation points PP is evaluated 
from the cell averages by a matrix-vector multiplication 
= Cfu.+r originating from ( |23[ ); (iv) the final re¬ 
construction is calculated by (30). For simple WENO 


computations, the only difference is that in steps (i) and 
(iii) the matrices are larger, because the Pi andpi(^9j (Jq 
not only depend on the averages of the neighbor cells, 
but the full polynomial given by TV -|- 1 coefficients per 
cell. 

In contrast to previous work, where no restriction al¬ 
gorithm (Sec. IIC) was present, we need to communi¬ 


cate more then just the two-dimensional bonndary lay¬ 
ers of every cell. Therefore, we introduced a new grid 
distribution method to reduce the communication be¬ 
tween different processors. The Cartesian grid consist¬ 
ing of n = UxTiyUz boxes is distributed on p processes 
in a way that communication between the processes is 
minimal. For this purpose, we perform a prime decom¬ 
position of p = PiP 2 ---Pi and set the number of grids per 
direction px = Pi,Py = P 2 ,Pz = P 3 initially. Let p^in = 


TABLE I: Numerical errors and convergence orders for the 
advection equation problem ( |46[ ) at t = 10 for different num¬ 
bers of grid patches and orders of DG polynomials N 
(CFL= 0.25, M = 1). 



DG 

DG-bWENO-7 

DG-bsimple WENO 

rix N 

Li error order Li error 

order L\ error 

order 

10 1 

1.53-10“^ - 

2.92-10“’ 

- 

1.51-10“’ 

- 

20 

3.63-10“^ 2.08 

1.40-10“’ 

1.05 

4.79-10“® 

1.65 

40 

5.34-10“® 2.76 

3.77-10“® 

1.90 

6.64-10“® 

2.85 

80 

7.29'10-‘ 2.87 

7.15-10“® 

2.39 

7.29-10“"* 

3.18 

160 

1.24-10“'‘ 2.54 

1.28-10“® 

2.48 

1.24-10“* 

2.54 

320 

2.84-10“® 2.12 

2.36-10“^ 

2.43 

2.84-10“® 

2.12 

10 3 

2.04-10“'‘ - 

4.84-10“® 

- 

2.12-10“* 

- 

20 

1.02-10“® 4.32 

1.51-10“® 

4.99 

1.02-10“® 

4.37 

40 

6.27-10“'^ 4.03 

4.99-10“® 

4.92 

6.36-10“'^ 

4.01 

80 

3.90-10“® 4.00 

9.71-10“’’ 

5.68 

3.98-10“® 

3.99 

160 

2.44-10“® 4.00 

1.60-10“® 

5.92 

2.53-10“® 

3.97 

320 

1.52-10“^° 4.00 

3.03-10“’° 

5.72 

1.62-10“’° 

3.96 

10 5 

7.99-10“’’ - 

9.95-10“® 

- 

5.75-10“° 

- 

20 

1.88-10“® 5.40 

1.11-10“® 

3.15 

1.27-10“'^ 

5.49 

40 

8.92-10“^° 4.39 

4.09-10“"’ 

4.76 

5.26-10“® 

4.59 

80 

5.42-10“!^ 4.04 

8.67-10“° 

5.56 

1.28-10“’° 

5.35 

160 

3.40-10“^® 3.99 

1.40-10“’’ 

5.95 

1.52-10“” 

3.07 

320 

6.48-10“^® 2.39 

7.66-10“’° 

7.51 

1.94-10“’° 

- 


UUTL(px,Py,Pz), we recalculate Pmin as Pmin Pmin • Pi- 
For further p^, we proceed in the same manner, so that 
each pj is always multiplied with the smallest of Px,Py,Pz- 
Finally, we subdivide the full box grid in px parts in x 
direction, in py parts in y direction and in pz parts in z 
direction. Each of these Px ■ Py ■ Pz parts is mapped to 
one MPI-process, which gives a simple box decomposi¬ 
tion, which is almost cubical. 


IV. SIMPLE TESTBEDS 


A. Advection equation 


As a first test for our new algorithms we consider the 
advection equation without a source {S = 0) 

dtu -I- dxU = 0 (45) 

for a Gaussian peak on the interval x £ [~1,1] 

V>(a;, 0 ) = 

(46) 

(we artificially add two peaks to gain smooth, periodic 
initial data) and a rectangular pulse (non-smooth initial 
data) 


ip{x,Q) 


1 if |a; — xoI < cr 
0 else 


(47) 
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The convergence rate in the first test case {A = 1, a = 
0.4) is influenced by several effects; see Table. |Tj 

For our choice of polynomials with order iV, we find 
convergence rates up to order + 1, as expected. How¬ 
ever, in an error regime beyond 10“we observe a fur¬ 
ther drop in the convergence rates, because of the grow¬ 
ing influence of truncation errors. Applying the stan¬ 
dard WENO reconstruction procedure leads to slightly 
different results. Convergence for small numbers of 
is slower, but finally shows convergence above N + 1-th 
order. This can be explained by the decreasing influ¬ 
ence of the WENO procedure for increasing rix- The 
cell indicator only marks the cells around the maxi¬ 
mum of the Gaussian peak as troubled, so the effective 
area, where the WENO reconstruction takes place, de¬ 
creases. Since the reconstruction has a strong smooth¬ 
ing effect, the numerical results significantly differ from 
the analytic solution for small Ux and tend to the pure 
DG solution for large Ux- Gomparing the two WENO- 
implementations, we observe that the simple WENO al¬ 
gorithm shows slower convergence, but while the stan¬ 
dard WENO is ~ 1 — 2 orders of magnitude less accurate 
than the pure DG-evolution, the simple WENO performs 
much better, showing roughly the same Li-errors as the 
pure DG-evolution. 

For the second case Eq. (471, which we just want to 
summarize briefly, we observe larger total errors than 
for the smooth problem discussed above. Again the pure 
DG-method errors are below the corresponding errors for 
the DG -I- standard WENO method. However, the dif¬ 
ference are at most a factor of 2. The simple WENO 
algorithm has comparable errors as the DG -I- standard 
WENO method. Independent of the scheme we observe 
first order convergence, which is consistent with the ex¬ 
pectation for a non-smooth problem containing disconti¬ 
nuities. 


B. Burgers equation 

The Burgers equation without source {S = 0) 


t 

0.00 0.05 0.10 0.15 0.20 0.25 0.30 



FIG. 2: Convergence rate during the evolution of a Gaussian 
wave packet for the Burgers equation (48l: As expected, the 
convergence rate is around A -|- 1 during the evolution of a 
smooth wave. At tahock ~ 0.233, a shock forms and the rates 
significantly drop down to first order convergence. When a 
standard WENO-7 reconstruction is used (crosses) the con¬ 
vergence rates are slightly higher than for the pure DG scheme 
(dots). The convergence rate is calculated from the errors of 
a. rix = 160 and a Ux = 320 run. 


up to tshock- Shortly before the shock formation at Ghock 
convergence start to drop for all N (gray shaded region). 
Employing a standard WENO algorithm convergence is 
slightly above the expected N + 1-th order. As discussed 
for the advection equation, this is related to the amount 
of troubled cells, which are reconstructed. For higher 
resolution a smaller percentage of cells is reconstructed 
and consequently a faster convergence is observed. After 
the shock formation the convergence order drops also for 
DG-hWENO to approximately first order convergence. 

In addition, we prepared the initial conditions 


dtu + udxU = 0 


(48) 


'0(a;, 0) = 0.5-|-sin(a;7r) (50) 


allows the formation of shocks from smooth initial data 
uq. After the time 


Ghock = - min 


duQ 

dx 


-1 


(49) 


shocks will appear during the evolution. We use this as 
a testbed for our code and evolve the initial Gaussian 


peak (461 with A = 1 and a = 0.2. For this initial condi¬ 
tions, a shock forms at Ghock ~ 0.23316. Similarly to our 
results for the advection equation we observe that the 
convergence rate is decreasing after tshock; cmp. Fig. 
The upper panels show snapshots of the field u at times 
t = 0.05,0.15,0.25. Without WENO-reconstruction (cir¬ 
cles) we observe the expected convergence order of A -|-1 


and check convergence at t = 0.5 /tt to compare with 
the results of [5], Tab. [H] summarizes the results. Be¬ 
cause of the smoothness of the solution, we observe for 
N = 1 polynomials second order convergence indepen¬ 
dent of the reconstruction method applied in the trou¬ 
bled cells. While the total Li-error for DG-l-WENO-5 is 
approximately a factor of 2-3 larger than the pure DG 
evolution, we see that the DG-l-simple WENO algorithm 
performs as good as pure DG. For A = 3 polynomi¬ 
als, we expect fourth order convergence, which we can 
verify with the pure DG and the DG-|-simple WENO 
setup. The DG-l-WENO-5 algorithm shows a higher con¬ 
vergence rate for low resolutions, which is again caused 
by the fact that larger number of cells decrease the inter¬ 
val where a reconstruction is performed. 
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TABLE II: Numerical errors and convergence orders for the 
Burgers equation problem (501 at t = ^ for different num¬ 
bers of grid patches rix and orders of DG polynomials N. 



DG 


DG -k WENO-7 DG -k simple WENO 

tlx N 

Li error 

order Li error 

order 

Li error 

order 

10 1 

5.34-10“^ 

- 

6.09-10“® 

- 

8.23-10“® 

- 

20 

1.45-10“^ 

1.87 

1.80-10“® 

1.75 

1.61-10“® 

2.34 

40 

4.29-10“^ 

1.76 

4.66-10“® 

1.94 

4.80-10“® 

1.75 

80 

1.24-10“® 

1.78 

1.29-10“® 

1.84 

1.24-10“® 

1.94 

160 

3.60-10“'‘ 

1.78 

3.69-10“’ 

1.80 

3.60-10“’ 

1.78 

320 

1.02-lO-'^ 

1.82 

1.03-10“’ 

1.83 

1.02-10“’ 

1.82 

10 3 

1.80-10“^ 


3.94-10“® 

- 

1.80-10“® 

- 

20 

9.80-10“® 

4.20 

1.50-10“’ 

4.71 

9.80-10“® 

4.20 

40 

6.36-10“® 

3.94 

6.72-10“° 

4.48 

6.36-10“° 

3.94 

80 

4.21-10-^ 

3.91 

4.22-10“’’ 

3.99 

4.21-10“’ 

3.91 

160 

2.71-10“® 

3.95 

2.71-10“® 

3.95 

2.71-10“® 

3.95 

320 

1.75-10“® 

3.95 

1.75-10“° 

3.95 

1.75-10“° 

3.95 

10 5 

3.58-10“® 

- 

5.86-10“® 

- 

3.58-10“® 

- 

20 

7.61-10“’’ 

5.55 

1.49-10“’ 

5.29 

7.60-10“’ 

5.55 

40 

1.61-10“® 

5.56 

1.33-10“° 

6.81 

1.62-10“® 

5.54 

80 

2.97-10“’° 

5.75 

1.16-10“® 

6.83 

2.98-10“’° 

5.77 

160 

5.47-10“’® 

5.76 

1.63-10“’° 

6.14 

5.47-10“’® 

5.76 

320 

1.15-10“’® 

5.56 

9.49-10“’® 

7.43 

1.15-10“’® 

5.57 


V. SPECIAL RELATIVISTIC 
HYDRODYNAMICS 

In the following section, we solve the GRHD conser¬ 
vation law 0 0 UHl without source terms and with 
q; = -0^ = 1 to consider special relativistic test cases, 
i.e. flat spacetimes. 


A. One-dimensional problems 


As a first test, we consider a smooth sine wave propa¬ 
gating with constant speed. The initial conditions are: 


p{x, t) =1 -I- 0.2 sin(27r(a; — Vxt)) 

Vx{x,t) = 0.2 

pix,t) =1 


(51) 


inside the periodic ID domain x G [—1)1] divided into 
Ux uniform grid patches. Viewing the Li errors and con¬ 
vergence rates (Tab. Ill), we find the convergence rate of 
the DG scheme to be V -f 1, when we use polynomials 
pGP^([-l,l]). 

Although we are dealing with a smooth problem a few 
cells around the maximum of the density p are marked 
as troubled. When we employ the standard WENO-5 
or WENO-Z reconstruction method, we observe at least 
one order of magnitude larger absolute errors as in the 
the pure DG-case for the employed resolutions. Contrary, 
the convergence order is artificially higher than for the 



FIG. 3: Evolution of the special relativistic shocktube ini¬ 
tial data ( |52[ | (density: red, velocity: blue, pressure: green): 
Numerical results using the standard WENO-3 (top, dots), 
standard WENO-5 (top, crosses), WENO-Z (bottom, dots), 
simple WENO (bottom, crosses) and the subcell evolution 
method (diamonds), compared to the analytical result (black 
line) at t — 0.4. For the troubled cell indication, we set 
M = 5. 


pure DG-method. For the simple WENO method, we 
obtain absolute errors compatible or identical with the 
scheme without reconstruction and obtain a convergence 
order oi N + 1 for an 7V-th order polynomial. In case of 
the subcell evolution, i.e. when we project the grid patch 
data on a finer subcell treating this with finite differenc¬ 
ing method, we observe similar convergence rates. The 
subcell evolution itself is performed with a fifth order 
accurate scheme [49j . which we verify with simulations 
using only subcells. 

As a second test focusing on the ability of our scheme 
to deal with discontinuities, we consider the shock tube 
problem with initial conditions 


{p,Vx,p){x,0) 


(10,0,13.33) ifx<0.5 
(1,0,10-^) if a: >0.5 


(52) 


on the domain x G [0,1]. The analytical solution for 
this problem in the context of SRHD is given by [50]. 
During our tests, we observe the troubled cell indica¬ 
tor to work reliable, since the grid patches which evolve 
the shock and the rarefraction wave are marked as trou¬ 
bled. All methods, the standard DG-WENO methods, 
the simple WENO approach as well as the subcell pro¬ 
jection method, are able to provide a stable evolution of 
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DG 


DG -b WENO-5 

DG + WENO-Z 

DG + simple WENO 

DG + subcells 

sub cells only 


N 

Li error 

order 

Li error 

order 

L\ error 

order 

Li error 

order 

Li error 

order 

Li error 

order 

10 

1 

1.22-10"® 

- 

8.46-10“® 

- 

8.62-10“® 

- 

1.80-10“® 

- 

2.85-10“® 

- 

3.22-10“® 

- 

20 


2.73-10“'* 

2.15 

2.80-10“® 

1.59 

2.69-10“® 

1.67 

1.86-10“® 

3.27 

1.95-10“® 

0.54 

1.00-10“* 

5.00 

40 


6.72-10“® 

2.02 

4.83-10“® 

2.53 

4.77-10“® 

2.49 

7.03-10“® 

4.72 

4.07-10“'* 

2.26 

3.14-10“° 

4.99 

80 


1.67-10“® 

2.00 

6.61-10“'* 

2.86 

6.43-10“'* 

2.89 

1.67-10“® 

2.06 

8.86-10“® 

2.20 

9.84-10“** 

4.99 

160 


4.18-10“® 

2.00 

9.64-10“® 

2.77 

8.73-10“® 

2.88 

4.18-10“® 

2.00 

2.02-10“® 

2.13 

3.08-10“*® 

4.99 

320 


1.04-10“® 

2.00 

1.40-10“® 

2.77 

1.44-10“® 

2.59 

1.04-10“® 

2.00 

4.31-10“® 

2.22 

1.14-10“*® 

4.75 

10 

3 

4.27-10“® 

- 

3.69-10“® 

- 

9.67-10“"* 

- 

4.33-10“® 

- 

1.58-10“® 

- 

1.08-10“* 

- 

20 


3.29-10“'^ 

3.70 

4.52-10“® 

6.35 

1.83-10“® 

5.72 

3.21-10“* 

3.75 

9.33-10“* 

4.08 

3.39-10“® 

4.99 

40 


1.79-10“® 

4.20 

7.37-10“'* 

5.93 

2.10-10“* 

6.44 

1.76-10“® 

4.18 

4.49-10“® 

4.37 

1.06-10“*° 

4.99 

80 


9.39-10“*® 

4.25 

1.09-10“® 

6.07 

3.39-10“° 

5.95 

9.50-10“*° 

4.21 

3.56-10“° 

3.65 

3.31-10“*® 

5.00 

160 


6.01-10“** 

3.96 

1.58-10“*° 

6.10 

1.31-10“*° 

4.69 

6.06-10“** 

3.96 

2.08-10“*° 

4.09 

1.07-10“*® 

4.94 

320 


3.80-10“*^ 

3.97 

6.96-10“*® 

4.51 

6.93-10“*® 

4.24 

3.84-10“*® 

3.97 

1.26-10“** 

4.04 

2.06-10“*'* 

2.37 

10 

5 

2.63-10“® 

- 

8.53-10“® 

- 

1.48-10“® 

- 

3.79-10“® 

- 

5.09-10“® 

- 

1.78-10“® 

- 

20 


3.86-10“** 

6.08 

2.67-10“'* 

4.99 

2.13-10“® 

6.11 

5.55-10“*° 

6.09 

1.36-10“° 

5.22 

5.57-10“*° 

4.99 

40 


6.13-10“*® 

5.97 

4.80-10“® 

5.79 

1.89-10“* 

6.81 

6.97-10“*® 

6.31 

1.04-10“** 

7.02 

1.74-10“** 

4.99 

80 


4.64-10“*^* 

3.72 

6.47-10“® 

6.21 

1.54-10“° 

6.93 

1.95-10“*® 

5.15 

3.93-10“*® 

4.73 

5.52-10“*® 

4.97 

160 


8.63-10“*'* 

- 

2.98-10“*° 

7.76 

1.30-10“** 

6.89 

9.70-10“*'* 

1.00 

1.20-10“*® 

1.70 

4.41-10“*"* 

3.64 

320 


1.80-10“*® 

- 

8.81-10“*® 

8.40 

6.47-10“*® 

4.32 

7.69-10“*® 

- 

2.48-10“*® 

- 

4.94-10“** 

- 


TABLE III: Numerical errors and convergence orders for problem (511 at t 
of DG polynomials N and several shock resolution methods. 


2 for different numbers of grid patches n^;, orders 


the shock tube problem, shown in Fig. 


B. Two-dimensional problems 

Generalizing our results to more complex two- 
dimensional wave setups, we perform two tests as pre¬ 
sented [23|: A shock-like test with the initial conditions 


(p,Vx,Vy,p)(x,0) = 


(0.03515,0,0,0.163) 

a X > Q,y > 

(0.1,0.7,0,1) 

a X < Q,y > 

(0.5,0,0,1) 

X < Q,y < 

(0.1,0,0.7,1) 

a X > Q,y < 


(53) 

with the initial conditions 

(0.5, 0.5,-0.5, 5.0) 

if a: > 0, y > 

(1,0.5,0.5,5.0) 

if a: < 0, y > 

(3.0,-0.5,0.5,5.0) 

if a < 0, y < 

(1.5,-0.5,-0.5,5.0) 

if a > 0, y < 


(54) 


(p,Vx,Vy,p)(x,0) = 


with (x,y) € [—1,1] X [—1,1]. During the evolution of 
both cases, all initial discontinuities are captured by the 
troubled cell indicator. We get the results as shown in 
Fig.i We tested in detail the standard WENO and the 
DG-bsubcell scheme, the figures show that the WENO- 
5 and DG-bsubcell evolution give qualitatively the same 
results. In case of the shocktube, Eq. (53)- left pan¬ 
els, less cells are marked troubled for the DG-|-subcell 


scheme. Furthermore, the DG-|-subcell method resolves 
steep gradients better than the standard WENO recon¬ 
struction. This becomes most dominant in a domain 
around x = y = —0.2. However, due to the larger com¬ 
putational expenses the DG-bsubcell scheme is a factor 
of ~ 2.4 times slower than the standard WENO method. 


The right panels of Fig. represent the vortex test, 
cmp. (54). As for the shocktube, both methods are able 
to resolve the structure properly. Again the DG-bsubcell 
method gives more accurate results, i.e. acting less dis¬ 
sipative keeping shock regions resolved, but also needs 
more computational resources and is ^ 3.2 times slower 
than the standard WENO implementation. 


VI. GENERAL RELATIVISTIC 
HYDRODYNAMCIS 


As the final test of our new implementation, we con¬ 
sider relativistic material in a curved spacetime back¬ 
ground and present results for a TOV-star in Cowling- 
approximation in ID, 2D, and 3D. Notice however, that 
the ID and 2D description is not identical to the 3D 
star. Being more specific, surfaces of constant densities 
correspond for the ID test to planes, for the 2D test 
to cylindrical shells, for the 3D test to spherical shells; 
cmp. discussion below. 
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FIG. 4: Special relativistic hydrodynamics simulations in 2D for the shocktube problem (531 (left) and the vortex problem 
(right) at t = 0.8, each evolved with the standard WENO-5 reconstruction (top) and the subcell evolution method (bottom) 
using n = 100 x 100 grid patches, polynomials of order N = 3, M = 5, CFL= 0.25. Density plot with contours, corresponding 
velocity field (arrows) and tronbled cells (shaded regions). 


A. Initial configuration 


Initial configurations for a single spherical symmet¬ 
ric neutron star are obtained by solving the TOV equa¬ 
tion [23112S]. The four-metric for a TOV star is given 
by 


ds^ = -e2<^dt2 -H 



-1 

dR^ + R^dn^. 


(55) 


To obtain m{R),(j){R), and the pressure p{R), the 
equations 

m + dirr^p 1 


dR 


— = {p{l + e)+p) 


R{R -2m) ^ ’ 


TOV 

(56) 


+ 

dcj) m + 4 :TtR^p 
cLR “ R{R-2mY 


(57) 

(58) 


are solved with an explicit fourth order Runge-Kutta al¬ 
gorithm. As starting values p{R = 0) = Pcentrab w(i? = 
0) = 0, and (f){R = 0) = 0 are specified and the system is 
closed by the polytropic EOS; Eq. §. Afterwards a coor¬ 
dinate transformation is performed to obtain the metric 
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FIG. 5: Density (blue), density error (green) and velocity (red) cell averages at t = 1000 for a ID TOV star using n = 100 
cells, polynomials of order N = 3, CFL = 0.25, fatm = 1-10“® and fthr = 100 evolved with DG and several shock resolution 
methods. The cells marked as troubled at f = 1000 are colored in gray. The result of the pure subcell run, which is equivalent 
to a finite difference simulation, is shown for comparison. 


in isotropic coordinates, which we use for the evolution, 
because a and ip‘^ can easily be obtained from this form. 
This solution describes the spacetime of a static, spheri¬ 
cally symmetric star. However, due to the discontinuity 
at the stars surface and truncation errors, the evolution 
is non-trivial. 


B. ID-TOV tests 


Before we are going to investigate the performance of 
our newly implemented algorithms in full 3D-simulations, 
we want to consider configurations similar to the TOV- 
star in just one dimension. We do not use spherical polar 
coordinates and stay in our Cartesian coordinate frame¬ 
work. Thus, all derivatives along the y- and z-direction 
are set to zero to achieve translation symmetry, i.e. dyf^ 
and in Eq. ([^ are zero and also all hrst derivatives 


in S Therefore, the obtained spacetime is different to 
a spherical symmetric TOV star. Nevertheless, it is still 
a valid testbed for our numerical scheme and with the 
restriction to a fixed spacetime background, the initial 
condition are in hydrodynamical equilibrium. 


Because of the smaller computational costs, we will 
discuss in detail ID-TOV results for all reconstruction 
algorithms, in particular we study WENO-3, WENO- 
5, WENO-Z, simple WENO reconstruction, as well as a 
DG-bsubcell and a pure subcell method for comparison. 
We have set in all our tests /atm = 10“® and /thr = 10^. 
Figurej^shows the density p (blue), the velocity (red), 
and the difference \p — Panail (green), where Panai refers 
to initial condition constructed according to Sec. VIA 


^ Notice that no second derivatives are present in Eq. 0 and 
that due to the restriction to Cowling approximation also first 
and second derivatives present in the metric field equations do 
not affect the simulation. 
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FIG. 6: Convergence test for a ID TOV star a.t t = 100 for three resolutions rihigh = 100, rimid = 50, niow = 25, polynomials of 
order N = 3, CFL = 0.25, fatm = 1-10“® and fthr = 100 evolved with DG and several shock resolution methods. 


All reconstruction algorithms lead to stable evolutions. 
In general we observe 3 regions of troubled cells, the 
left star surface, the maximum of the density, and the 
right star surface. During the evolution some troubled 
cells are activated or deactivated, which explains why for 
WENO-Z reconstruction at the presented time t = 1000 
the surfaces are not marked as troubled. 

We observe that WENO-3, WENO-5, WENO-Z per¬ 
form worst, i.e. large velocities are present at the stars’ 
surface and |p—Panail is larger as for the other reconstruc¬ 
tion mechanisms (notice the different y-scales for and 
|p — PanaiD- The best results are obtained with the simple 
WENO and DG-|-subcell methods. The total Li errors 
of p for the given setup are 6.0T0“^ for simple WENO 
and 4.6-10“^ for DG-|-subcell method. The pure subcell 
evolution performs as good as the DG-|-subcell method. 

The advantage of the simple WENO and subcell meth¬ 
ods can be understood by considering the effectively 
higher resolution compared to the other schemes. In the 
standard WENO case, only the cell averages are used 
for the componentwise reconstruction, therefore the ef¬ 
fective resolution drops depending on the employed poly¬ 
nomial order. In contrast, the simple WENO approach 
uses the full information of the polynomial inside the 
cell and additionally uses only three cells for the recon¬ 


struction, thus no significant performance loss is obtained 
and the simple WENO reconstruction is a factor 1.57 
slower than the standard WENO-3 approach (a factor 
1.40 slower than the standard WENO-5 approach). Fi¬ 
nally, in the DG -I- subcell method points are added in 
problematic regions. Because of the additional compu¬ 
tational effort due to the projection between top- and 
subcells and the larger number of points in the troubled 
cells, the algorithm is a factor of ~ 1.67 slower than the 
standard WENO method. Although not noticeable for 
ID setups, we encounter for higher dimensional setups a 
significantly larger amount of memory, i.e. a ^ 2.7 times 
higher memory load for 2D runs (~ 4.8 times higher for 
3D) when subcells are activated compared to standard 
WENO-3 simulations. Nevertheless the DG-|-subcell ap¬ 
proach seems to be a valid choice for further development, 
while it allows (i) to reuse well-tested ED schemes in trou¬ 
bled regions, (ii) give the most accurate results due to an 
effectively higher resolution in troubled regions, (iii) al¬ 
lows a speed up compared to the usually employed ED 
codes, because of a more effective DG method in large 
parts of the numerical domain. 

In Fig. we present a pointwise convergence test for 
all methods. We compare evolutions with 25, 50,100 cells 
and use polynomials of order 3. The difference between 
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FIG. 7: Convergence order for a ID TOV star during evolu¬ 
tion t £ [0,100] for two resolutions Whigh = 100, niow = 50, 
polynomials of order N = 3, CFL = 0.25, fatm = I-IO"® and 
fthr = 100 evolved with DG and three shock resolution meth¬ 
ods: standard WENO-3 (top), simple WENO (middle) and 
DG-fsubcell (bottom). 


the low and medium resolution is shown blue, while the 
difference between the medium and high resolution is 
shown red. We rescale the difference of the medium and 
high resolution according to the expected convergence 
order, i.e. 3rd order for WENO-3 and 4th order for the 
other schemes. We observe that in all cases we obtain 
roughly the expected convergence order. Furthermore in 
the logarithmic plots is clearly visible that for some se¬ 
tups the outer regions of the star is smeared out. In case 
of standard WENO algorithms larger stencils (WENO-5 
and WENO-Z) lead to a numerical solution where the 
outer star layers are not fixed and no sharp surface is 
visible, this improves for the WENO-3 reconstruction. 
Contrary, the simple WENO and DG-|-subcell method 
keep the surface of the star fixed. In all runs higher res¬ 
olution improves the results and less material is leaving 
the star. 

The simplicity of the ID-TOV star allows us to con¬ 
sider setups with higher resolution than achievable in the 
corresponding 2D and 3D tests and a more detailed anal¬ 
ysis becomes possible. A recurring question is how re¬ 


gions with low order convergence (because of low differen¬ 
tiability of the solution) affect regions where the solution 
is smooth. Specifically, do the regions of low order remain 
localized, or if not, how quickly does the loss of conver¬ 
gence spread through the entire domain? See for example 
[51) . where the wave equation with discontinuous initial 
data is studied, for which analytic results are available 
in [52] predicting the growth of the non-convergent area 
with, e.g., the square-root of time, ^ '/i. 

Figure shows the convergence order during the first 
stages of the evolution for 50 and 100 cell setup. Pre¬ 
sented are the WENO-3 (top panel), simple WENO (mid¬ 
dle panel), and the subcell (bottom panel) evolution. For 
all panels, we observe that inside the star, where also 
cells are marked as troubled, the WENO-3 method shows 
~3rd order convergence and the simple WENO method 
a convergence order above 4. Furthermore, while for 
WENO-3 the error seems to corrupt the convergence in 
the entire star it seems to be localized for simple WENO 
for the entire simulation. For the subcell evolution we 
observe that the convergence order at the stars’ center 
lies between second and third order, which is consistent 
with the employed flux methods implemented in the FD 
subcells [15]. Artificially setting the center cells non- 
troubled cures this problem and leads locally to higher 
order convergence. However, it has no influence on the 
global convergence order. More problematic, a large error 
is traveling inwards from the outer surface for all simu¬ 
lations, which leads to a lower convergence order for all 
setups. It is important to notice that this effect is not 
related to the movement of troubled cells. The region 
of troubled cells stays relatively fixed at the stars’ sur¬ 
face. However, the flux across the cell surfaces seems to 
contain lower order components. Regarding this fact, it 
is debatable whether one can obtain high order conver¬ 
gence in more general setups, e.g. dynamical spacetimes 
and moving objects. 


C. 2D TOV star 

Considering the results of the previous section, the sim¬ 
ple WENO and DG-|-subcell schemes seem to be prefer¬ 
able. However, we found, that the simple WENO method 
performs worse in higher dimensional problems as in the 
ID case. Compared to the standard WENO reconstruc¬ 
tion, where a smoother polynomial from several cell av¬ 
erages is constructed, the simple WENO methodology 
allows steeper gradients and has weaker smoothing influ¬ 
ence. For runs of higher dimensional problems we ob¬ 
serve this smoothing to be crucial for the stability of the 
evolution. Furthermore, the simple WENO computation 
underlies a significant slowdown in d > I Dimensions, 
because the evaluation of the smoothness indicators is a 
quadratic form of all {N + 1)‘^ coefficients. This is the 
reason why the standard WENO-3 scheme, which turns 
out to allow stable evolutions, is used instead. 

We investigated the convergence of the two schemes. 
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FIG. 8: Density Li error for a 2D TOV star at t = 500, 
t — 1000 for six resolutions n = 25, 50, 75,100,150, 200, poly¬ 
nomials of order N = 3, CFL = 0.25 and fthr = 100 evolved 
with DG -I- standard WENO-3 and the DG -I- subcell evo¬ 
lution method. For WENO-3 we set fatm = I-IO”®, for the 
DG -I- subcell we set fatm = I-IO”®. The dashed black lines 
correspond to second / third order convergence. 
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simple WENO and DG-|-subcell, for the 2D TOV star^ 
regarding the density Li error. As shown in Fig. we 
observe a convergence order of ~ 2 for the subcell scheme, 
which indicates that the evolution error originating in the 
subcells spreads over the grid and leads to a lower order 
of convergence. In comparison, the standard WENO-3 
scheme converges in third order for coarse grids. The 
subcell method causes a much higher memory load and 
longer calculation times (see Sec. VIB) because of the 
higher number of grid points in each direction. For this 
reason, we decided to use the standard WENO-3 method 
for the 3D simulation of a TOV star. 


D. 3D TOV star 

Considering a 3D TOV star, we are able to provide 
a stable simulation with a DG -I- standard WENO-3 
method. Although we show our results up to t = 1000, 
there is no evidence of any instabilities for longer runs. 
We are considering two numerical setups at resolutions 
25 X 25 X 25 and 50 x 50 x 50 and the initial configu¬ 
ration as a reference solution. After a short transition, 
the numerical simulations reach an almost steady struc¬ 
ture as shown in Fig.|^ The density prohles along the x- 
and y-axis are shown as red and blue lines. The differ¬ 
ence between the densities for z = 0 is presented as the 
gray shaded region. On the bottom panel, we present the 


2 As for the ID test, we employ Cartesian coordinates and due to 
the restriction to a fixed spacetime background also our 2D-TOV 
example is in hydrodynamical equilibrium. 



FIG. 9: Pointwise convergence order in the z = 0 plane for the 
density of a 3D TOV star at t = 500 using two resolutions 
tihigh = 50, niow = 25, polynomials of order A = 3, CFL 
= 0.25, fatm = I-IO”® and fthr = 100 evolved with DG and 
standard WENO-3. The two density solutions for n — 50 
(red) and n = 25 (blue) are shown on the axes x = z = 0 
and y = z = 0. Their difference in the z-plane is shown in 
gray, the corresponding zero-crossing is indicated by the gray 
contour line. 
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FIG. 10: Density Li error for a 3D TOV star at t = 500, 
t = 1000 for four resolutions n = 20, 32, 40, 50, polynomials 
of order A = 3, CFL = 0.25, fatm = I-IO”® and fthr = 100 
evolved with DG and standard WENO-3. The dashed black 
lines correspond to third order convergence. 


computed convergence order. In large areas of the star 
second to fourth order convergence is present and even 
higher convergence in its center and outside areas near 
the surface. The latter can be explained by the failure of 
the coarse grid setup to keep the density on atmosphere 
level outside the star, whereas the fine grid setup does. 
The narrow band of low convergence (colored blue) inside 
the star can be explained as follows: The finer resolved 
solution stays closer to the density maximum in the star 
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center and zero at the star’s surface, the opposite holds 
for the coarse resolution. Thus, the differences tend to 
zero, see solid black line, and the convergence drops lo¬ 
cally. 

As a global measurement of the convergence order, we 
present the Li-norm for the 3D TOV star in Fig. 10 for 
four different resolutions. Similar to the 2D test case, we 
observe an almost third order convergence (black dashed 
line) for the standard WENO-3 algorithm with third or¬ 
der polynomials. Using a htting function of the form 
A ■ n~^ for the Li-error, the obtained convergence order 
is & = 2.75 for t = 500 and b = 2.88 for t = 1000 and 
thus close to the theoretical expected value. 


VII. CONCLUSION 

In this work, we presented new algorithms imple¬ 
mented in the existing bamps code: a DG, a WENO- 
DG, and a mixed ED -|- DG-algorithm combined with 
standard WENO [8l|45] and a simple (compact) WENO 
scheme [55]. We tested all algorithms and reconstruction 
methods with a number of tests starting with the advec- 
tion and Burgers equation, the main results being exam¬ 
ples for special and general relativistic hydrodynamics. 
In almost all cases, we were able to obtain the expected 
convergence order for smooth solutions and also found a 
proper shock treatment in case of jumps and discontinu¬ 
ities. 

Our main result was the simulation of a single TOV- 
star, which we modeled in the Gowling approximation, 
i.e. for static geometric variables. In fact, while it has not 
been attempted yet to apply the existing DG methods for 
the vacuum Einstein equations 1101121] to 3D GRHD, we 
have demonstrated recently that the pseudospectral mul¬ 
tipatch method of bamps |14j works well for demanding 
3D vacuum spacetimes (highly non-linear gravitational 
waves that collapse to a black hole). The pseudospectral 
penalty method developed in m can be viewed as a spe¬ 
cial case of a DG method for the full Einstein equations 
in a non-flux form. Furthermore, the present work on 
GRHD and the wave-collapse simulations are compatible 
in the type of variables and equations they use, and can 
run with the same spectral element grid and polynomial 
basis functions (Ghebyshev or Legendre Gauss-Lobatto 
grids). Therefore, we do not expect any immediate ob¬ 
stacle to combine the existing geometry code with the 
new GRHD methods. 

One simplifying restriction in our implementation was 
the usage of a simple troubled cell indicator. We intend 
to study the influence of different and more sophisticated 
troubled cell indicators. While our simple setup allowed 
an easy implementation and stable evolutions, it also 
marked maxima as troubled, which should be avoided 
in the future application of the code. 

Keeping the number of employed cells fixed, we found 
that for ID problems the subcell and simple WENO al¬ 
gorithms were the most accurate ones. This can be easily 


understood, since the standard WENO method is based 
only on cell averages for the reconstruction. In contrast, 
the simple WENO method uses the knowledge of the en¬ 
tire polynomial, and the subcell methods resolved trou¬ 
bled cells with effectively 2N -|- 1-times more points. In 
our examples, the simple WENO and subcell methods 
have some drawback for higher dimensions. Both meth¬ 
ods come with a significant overhead, and it is planned to 
investigate more efficient implementations in the future. 
More importantly, the direct application of the simple 
WENO reconstruction in 2D led to unexpected instabili¬ 
ties, for example in the computation of the primitive from 
the conservative variables, which one should be able to 
avoid. This issue certainly deserves further study since 
the simple WENO method is very promising based on 
the ID results. 

Due to the large computational cost of the subcell 
method, we only employed the standard WENO recon¬ 
struction in 3D and investigated the subcell method in 
2D. In our 2D examples, the subcell method turned out 
to be (as expected) approximately second order. For the 
standard WENO-DG method we found 3rd order con¬ 
vergence for low and second order convergence for high 
resolutions. The observed third order convergence is con¬ 
sistent with our results in the full 3D simulation. How¬ 
ever, for high number of cells we noticed that a higher 
than second order convergence in the matter variables 
seems to be hard to obtain, since (i) the computation of 
the Li-norm of the error emphasizes inaccurate, prob¬ 
lematic regions, and (ii) errors propagate from the sur¬ 
face of the star through the neutron star and “corrupts” 
the order of convergence. Although this fact can be seen 
as a setback, DG methods allow a better parallelization 
and refinement strategy than fixed FD codes and are one 
of the most promising methods to take into account for 
future GRHD-code developments. 

Our work can be seen as a first step towards a com¬ 
plete 3D-DG implementation for GRHD, since it employs 
DG-methods for GRHD problems beyond the limitation 
of spherical symmetry as in previous work. It is planned 
to further develop the numerical techniques by consider¬ 
ing adaptive mesh refinement, and to extend the physics 
to full general relativity (beyond the Cowling approxi¬ 
mation), which together will allow numerical simulations 
of astrophysical systems consisting of single and binary 
neutron stars. 
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Appendix A: DG method for scalar conservation 
laws in ID 


Following HIS], we summarize some aspects of the DG 
method that already arise in the non-linear, scalar, one¬ 
dimensional case. We add some details relevant to the 
present work concerning implementation issues and the 
equivalence of the once and twice integrated form of the 
equations. One of our goals is the combination of the 
DG method for relativistic matter with the pseudospec- 
tral penalty method of [T3] for the geometry, which is 
not using a flux conservative form, but is close to the 
strong formulation of the DG method given below in 
(A7). Therefore we examine the question how the dis¬ 


cretized equations for the weak and strong form are re¬ 
lated. 


1. Derivation of the discretized equations 

Gonsider 


dtu + d^f{u) = 0 (Al) 

for a function w(t, x) and a flux function f{u{t, x)) on the 
interval I = [—1,1]. Given a space of test functions on 
/, we obtain the weak form of the conservation law by 
integration. For a test function v(x), 

(v, dtu) + (v, dxf) = 0. (A2) 


Later we assume that the scalar product of two functions 
is (/i,/ 2 ) = fi{x)f 2 {x)dx, i.e. we assume the trivial 
measure which is the natural weight for the polynomial 
basis of Legendre polynomials. 

Part of the DG method is a special treatment of the 
flux at the boundaries of /, where we replace the flux / by 
a non-unique choice of a numerical flux /*. Integrating 
(A 2 ) by parts in space we arrive at two versions of the 
conservation law. 


{v,dtu) - = -[vf*], (A3) 

{v,dtu) + {v,da:f) = [v{f-f*)], (A4) 


which £i{xj) = 6 ij. We choose the Xi to be the colloca¬ 
tion points for Legendre-Gauss (LG) or Legendre-Gauss- 
Lobatto (LGL) integration. This allows the approxima¬ 
tion of u{x) by a nodal expansion that has the interpo¬ 
lation property u{xi) = Ui, Imu[x) = Yh=.q udi{x). 

When the meaning is clear from context, we write u 
instead of ImU- A key feature of the nodal expansion 
is that it works equally well for linear and non-linear 
functions, in particular 

N N 

u{x) = Uiti[x), f(u{x)) = ^ Mi{x), (A5) 
2=0 2=0 


where /* = f{u^) = fiu{xi)). 

The nodal approximation with iV-th-order polynomials 
leads to d iscretized versions of the conservation laws (A3) 
and (A4). Ghoose test functions v{x) = ii{x), and insert 
(A5) to obtain 


Mdtu-S^f = -[if*], (A 6 ) 

Mdtu + Sf = m-m (A7) 


where we have introduced the mass matrix M and stiff¬ 
ness matrix S, 


= = (A 8 ) 

We use matrix notation and a summation convention, 
e.g. Sf = Sij fj = J2f=p Sij fj ■ 

An important point is that in general the mass matrix 
is not diagonal, that is, the characteristic Lagrange poly¬ 
nomials are not necessarily orthogonal. Specifically, for 
LGL the mass matrix is not diagonal, while for LG it is 
diagonal. However, for both LGL and LG the matrix is 
symmetric and invertible. The stiffness matrix is directly 
related to the derivative matrix. 


Dij = da;£j{xt), (A9) 

which approximates the pseudospectral derivative at the 
nodes by {dxu){xi) = DijUj. We have [1] 

S = MD, M-^S = D, M-^S^ = 

(AlO) 

Given M, D, and a prescription for /*, we solve the 
explicit time-integration problem based on ( |A 6 [ ) or 

dtu-{M-^D'^M)f = -M-^[£f], (All) 
dtu + Df = M-i[^(/-/*)], (A12) 


where [g] = ( 7 ( 1 ) — g(—1) for any function g{x). Eqn. 
(A3) is obtained by integrating b y pa rts and replacing 
/ by /* at the boundary. Eqn. (A4) is obtained from 
(A3) by integrating by parts once more but leaving the 
resulting boundary term unchanged. We refer to (A3) 
and (A4) as the weak and strong form, respectively, or 
as the once and twice integrated flux equation to avoid 
confusion with the original “strong” form, (ig. 

The nodal DG spectral element method is based on a 
choice of A' -I- 1 distinct nodes Xi € /. Such nodes define 
the unique A-th-order Lagrange polynomials ii(x), for 


for the descretized, time-dependent function values Ui (t). 

The method generalizes immediately to a partition of 
any interval [a, 6 ] € K into several elements Ij with an 
appropriate mapping of the coordinates and with a cou¬ 
pling of neighboring elements through /*. 


2. Implementation issues 

Let us comment on some implementation issues, specif¬ 
ically for the LGL method. The nodes Xi and the LGL 
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integration weights Wi are obtained from the Legendre 
polynomials, for which simple but stable and accurate 
algorithms are available, e.g. [5]. The nodes Xi are the 
N — 1 roots of dxPN{x) combined with the endpoints of 
the interval, —1 and 1, for a total of iV + 1 nodes. The 
integration weights are Wi = 2/(TV(TV + l)PAr(xi)^). Var¬ 
ious other quantities are determined without further ref¬ 
erence to the Legendre polynomials by general formulas 
for Lagrange interpolation. The weights for barycentric 
interpolation are Ci = ^/{xi — Xj). The deriva¬ 

tive matrix is 


Dij = d^£j{xi) = 
N 

Pii — ^ ^ ■ 

T=0 


* h 


(A13) 

(A14) 


from the expansion of the Legendre polynomials in the 
Lagrange basis. Computing the difference to the diagonal 
approximation we find for LGL 


Ml, 


-1 


= —<5,: 


N+l 


PN{Xi)PN{Xj). 


(A19) 


Alternatively, note that M can also be computed directly 
as the analytic integral either by term by term 

integration after expanding the product ii{x)£j{x), or by 
exact Gauss integration on a secondary grid with TV -|- 2 
points. (In experiments, TV -|- 3 gives somewhat more ac¬ 
curate results.) However, we still have to find the inverse 
of M numerically. For large TV, (A19) may be preferred. 


3. Equivalence of once and twice integrated forms 


The equation for the diagonal term ensures that the 
numerical derivative of a constant like = 1 is zero 
|53| . Since the endpoints are included among the nodes, 
xq = —1 and xat = 1, 

[^igW = iiil)g{l) - £i{-l)g{-l) = SmgN - Siogo- 

(A15) 

There are several ways to compute the mass matrix 
Mij = {£i,£j)- One option is to perform the integra¬ 
tion numerically according to the Gauss formula associ¬ 
ated with the nodes, which approximates the integral of 
a function g{x) using the integration weights Wi, 

„i N 

/ g{x)dx o±'^Wig{xi). (A16) 

i=0 

This integration is exact if g{x) is a polynomial of degree 
up to 2TV -I- 1 for LG and up to 2TV — 1 for LGL. Since 
the integrand £i£j for the mass matrix is of degree 2TV, 
for LG the numerical integral is exact. 


Mij — {£i^£jl — {£ii £jl N — 


(A17) 


where {f,g)N = '^i'Wifigi denotes the numerical scalar 
product. However, for LGL we only obtain the approxi¬ 
mation 


Mij = {£i,£j) ~ {£i,£j)N = wA 


(A18) 


It turns out that this approximation, also called mass 
lumping, is equivalent to a certain filter that strongly af¬ 
fects the highest mode in the Legendre basis and that can 
reduce the effective order of the approximation |48j . In 
the context of spectral element methods of comparatively 
high order, say TV = 10, a pprox imating M for LGL by 
the diagonal matrix as in (AI8) is considered standard 
in 1^. However, for orders around TV = 2,3,4, it is of¬ 
ten preferable to evaluate Mij = {£i, £j) for LGL without 
approximation [4]. For example milH], TW-i = VV^, 
where V is the generalized Vandermonde matrix for the 
normalized Legendre polynomials. This relation follows 


For the continuum problem, we perform the integra¬ 
tion by parts 


{v,dxf) = [vf] - if,dxv). 


(A20) 


Under specific but quite general conditions the dis¬ 
cretized equations satisfy the corresponding summation 
by parts property exactly. In this case the once and twice 
integrated DG methods are numerically identical. There 
may be round-off errors, but there are no systematic er¬ 
rors that only converge away with increasing TV. This is 
fully explained in [IHl [Ml US] . 

Given the present setup, it is straightforward to show 
algebraic equivalence of (A6| and (A7). The difference 


between those two equations is 

Sf = [£f]-S^f, 
s = \££]-S^, 


(A21) 

(A22) 


for all fi, and independently of the choice of f* or the 
computation of M. In the transition to (A22) we use 


that / is approximated by an TV-th-order polynomial, 
(A5). By definition of Sij , 




— l£ii£dx£jl 


= [£r£j\-{£j,dx£i) 


(A23) 


so the summation by parts property (A22) does indeed 


hold. Summation by parts is exact for LG and LGL even 
if Sij is defined by numerical integration because 


S^, = {£t,dx£j) = {£t,dx£j)N 


(A24) 


since £idx£j is a polynomial of degree 2TV — 1. 

It is instructive to make the summation by parts for¬ 
mula for the LGL method more explicit. From (A241 and 


(A16) 


Sij — ^ ^^ '^k£i{xi£)dx£j{x}£^ — WiDij. (A25) 
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(Incidentally, this means th at Sg = MikD^j = WiDij for 
both LG and LGL). Hence (A22) becomes 




(A26) 


In other words, the summation by parts rule implies for 
LGL points a relation between the integration weights Wi 
and the barycentric interpolation weights c^, 


We now restrict ourselves to the LGL case. A priori it is 
not clear how a simple rescaling and a transpose of the 
derivative matrix leads to the term = {SiN — Sio)6ij, 
which is a diagonal matrix with non-vanishing entries 
only in two of the corners. For i = j, 


Di 

QiWiDi 


dxkixi) = T^{6iN - Sm), (A27) 

2^Wi 

6 ^N-SM = [^^^, (A28) 


so (A261 is satisfied on the diagonal. In particular, we see 
how the boundary terms come about. For i ^ j, (A26) 
becomes 


Wi 




from which we obtain with (A13l that 


Cj j 

^ i 


(A29) 


(A30) 




(A31) 


for some constant C. Surprisingly, the explicit relation 
between wfand was only found recently, see [SB] 
on such relations for Jacobi polynomials. For our case, 


cfGL = (A32) 


where Cat is an explicitly known constant that depends 
on the number of points. In summary, for LGL (or anal¬ 
ogously for LG), we can start from the general result 
on summation by parts and arrive at a partial proof of 
(A321, or we can start from relations like (A32) and prove 
summation by parts without directly using partial inte¬ 
gration in the continuum. 
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